import sys
sys.path.append("..")
import utils.common as common

import gdal
import osr
import gdalnumeric
import numpy as np
from numpy import linalg as LA
import cv2
import math
import os

def get_raster_info(raster_name):
	"""
	Args:

	Returns:

	Raises:
		RuntimeError: 
	"""
	if isinstance(raster_name,str):
		if not os.path.exists(raster_name):
			raise RuntimeError("Can not found %s"%raster_name)

		raster = gdal.Open(raster_name)
	else:
		raster = raster_name
	return __get_dataset_info___(raster)

def open(raster_name):
	"""
	Args:

	Returns:

	Raises:
		RuntimeError: 
	"""
	if not os.path.exists(raster_name):
		raise RuntimeError("Can not found %s"%raster_name)

	raster = gdal.Open(raster_name)

	return __get_dataset_info___(raster),raster.GetRasterBand(1).ReadAsArray()

def save(out_filename, data, fill_value, geoTrans, projection):
	"""
	Args:

	Returns:

	Raises:
		RuntimeError: 
	"""
	print("wirte to the %s" % (out_filename))

	# #covert the numpy type code to gdal type code
	# gdal_type = gdalnumeric.NumericTypeCodeToGDALTypeCode(data.dtype.type)

	# # check the return gdal type code
	# if type(gdal_type) != np.int:
	#     if gdal_type.startswith('gdalconst_GDT_') == False:
	#         gdal_type = eval('gdalconst.GDT_' + gdal_type)
	common.is_exists(os.path.split(out_filename)[0])

	if "int8" in data.dtype.name:
		gdal_type = gdal.GDT_Byte
	elif "int" in data.dtype.name:
		gdal_type = gdal.GDT_UInt16
	else:
		gdal_type = gdal.GDT_Float32

	# the data to save must be 2-dimension data
	if len(data.shape) != 2:
		raise ValueError('the out data must be 2-dimension numpy array!!')
		sys.exit(1)

	# replace the nan data
	data[np.isnan(data)] = fill_value

	# construct the out dataset
	out_ds = gdal.GetDriverByName('GTiff').Create(out_filename, data.shape[1],
												  data.shape[0], 1, gdal_type)
	out_ds.SetGeoTransform(geoTrans)

	out_ds.SetProjection(projection.ExportToWkt())

	# wirte the data to out dataset
	out_ds.GetRasterBand(1).WriteArray(data)
	out_ds.GetRasterBand(1).SetNoDataValue(fill_value)

	out_ds.FlushCache()

def detecte_isolate_pixel(img,window_size=3):
	"""
	Args:

	Returns:

	Raises:
		RuntimeError: 
	"""
	if not img.dtype == np.float:
		img = img.astype(np.float)

	im_fk3 = cv2.filter2D(img,-1,np.ones((window_size,window_size),np.float),borderType = 0)
	im_fk5 = cv2.filter2D(img,-1,np.ones((window_size+2,window_size+2),np.float),borderType = 0)

	## threshold
	return (im_fk3 > math.floor(window_size/2)) * (im_fk5 > math.floor((window_size+2)/2))

def latlont2imgXY(latlon,geo_trans):
	"""
	Args: latlon: [lon,lat]

	Returns: [imgX, imgY]

	Raises:
		RuntimeError: 
	"""
	return LA.solve([[geo_trans[1],geo_trans[2]],[geo_trans[4],geo_trans[5]]],[latlon[0]-geo_trans[0],latlon[1]-geo_trans[3]])

def imgXY2latlon(imXY,geo_trans):
	"""
	Args: [imgX,imgY]

	Returns: [lon,lat]

	Raises:
		RuntimeError: 
	"""
	return [geo_trans[0] + imXY[0]*geo_trans[1]+imXY[1]*geo_trans[2],geo_trans[3] + imXY[0]*geo_trans[4]+imXY[1]*geo_trans[5]]

def resample(src,ref,interp_method):
	"""
	Args: [imgX,imgY]

	Returns: [lon,lat]

	Raises:
		RuntimeError: 
	"""
	if isinstance(src,str):
		src_dst = gdal.Open(src)
	else:
		src_dst = src

	if isinstance(ref,str):
		ref_dst = gdal.Open(ref)
	else:
		ref_dst = ref

	out_dst_driver = gdal.GetDriverByName("MEM")
	out_dst=out_dst_driver.CreateCopy("",ref_dst)

	wrp_opts = gdal.WarpOptions(resampleAlg=__get_interp_method__(interp_method), srcNodata=src_dst.GetRasterBand(1).GetNoDataValue(), dstNodata=src_dst.GetRasterBand(1).GetNoDataValue(), multithread=True)
	# wrp_opts = gdal.WarpOptions(resampleAlg=__get_interp_method__(interp_method), srcNodata=src_dst.GetRasterBand(1).GetNoDataValue(), dstNodata=src_dst.GetRasterBand(1).GetNoDataValue(), multithread=True, callback=gdal.TermProgress)
	gdal.Warp(out_dst,src_dst,options = wrp_opts)

	return out_dst


def __get_interp_method__(interp_method):
	"""
	Args: [imgX,imgY]

	Returns: [lon,lat]

	Raises:
		RuntimeError: 
	"""
	interpolation_method = gdal.GRA_NearestNeighbour
	if interp_method == 'bilinear':
		interpolation_method = gdal.GRA_Bilinear 
	elif interp_method == 'lanczos':
		interpolation_method = gdal.GRA_Lanczos 
	elif interp_method == 'cubicspline':
		interpolation_method = gdal.GRA_CubicSpline 
	elif interp_method == 'nearestneighbour':
		interpolation_method = gdal.GRA_NearestNeighbour 
	elif interp_method == 'cubic':
		interpolation_method = gdal.GRA_Cubic
	elif interp_method == 'average':
		interpolation_method = gdal.GRA_Average
	elif interp_method == 'mode':
		interpolation_method = gdal.GRA_Mode
	else:
		raise Exception("Interpolation method was not recognised or known.")

	return interpolation_method

def __get_dataset_info___(raster_layer):
	"""
	Args:

	Returns:

	Raises:
		RuntimeError: 
	"""
	projection = osr.SpatialReference()
	result = projection.ImportFromWkt(raster_layer.GetProjectionRef())
	# if the dataset has not contian the SpatialReference
	# specific the SpatialReference as WGS84

	if result != 0:
		projection.SetWellKnownGeogCS('WGS84')

	return [
		raster_layer.GetRasterBand(1).GetNoDataValue(),
		raster_layer.GetGeoTransform(), projection
	]
